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Presented in this paper is the description of a Markov chain Monte Carlo (MCMC) routine for 
conducting coherent parameter estimation for interferometric gravitational wave observations of an 
inspiral of binary compact objects using multiple detectors. Data from several interferometers are 
processed, and all nine parameters (ignoring spin) associated with the binary system are inferred, 
including the distance to the source, the masses, and the location on the sky. The data is matched 
with time-domain inspiral templates that are 2.5 post-Newtonian (PN) in phase and 2.0 PN in 
amplitude. We designed and tuned an MCMC sampler so that it is able to efficiently find the poste- 
rior mode(s) in the parameter space and perform the stochastic integration necessary for inference 
within a Bayesian framework. Our routine could be implemented as part of an inspiral detection 
pipeline for a world wide network of detectors. Examples are given for simulated signals and data 
as seen by the LIGO and Virgo detectors operating at their design sensitivity. 

PACS numbers: 02.70.Uu, 04.80.Nn, 05.45.Tp, 07.05. Kf, 97.80.-d. 
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I. INTRODUCTION 

The era of gravitational wave astronomy is now close 
upon us as numerous interferometric detectors are op- 
erating. The Laser Interferometer Gravitational Wave 
Observatory (LIGO) 0, 0, EJ 

has now reached its target 
sensitivity, and there is the hope that a detection could 
come at any time [||. Around the globe a world-wide 
network of detectors is coming on-line; Virgo in Italy 




, GEO in Germany [8, 9], and TAMA in Japan 
11[ are operating alongside LIGO in the quest for 
gravity wave detection. These ground based laser inter- 
ferometers are sensitive to gravitational radiation in the 
frequency band from 40 Hz up to 8 kHz. 

Coalescing binaries containing neutron stars or black 
holes promise to be one of the cleanest and most prob- 
able sources of detectable radiation Observation of 
inspiral events could provide important information on 
the structure of neutron stars [lj, E3 • Even cosmologi- 
cal information can be extracted from the observation of 
inspiral events 0, El El El El- 

The characteristics 
of radiation in the post-Newtonian regime will provide 
insight into highly non-linear general relativistic effects, 
such as the observation of the formation of a Kerr black 
hole as the binary system decays [H, HO, Hl[ . The LIGO 
Scientific Collaboration (LSC) has been actively search- 
ing for binary inspiral events [22|, [23| , as well as conduct- 
ing searches in coincidence with TAMA 24 [. Using the 
data from LIGO's S2 run, it was possible to set an upper 
limit on the neutron star coalescence rate of less than 
50 per year per Milky Way equivalent galaxy 23]. The 
LSC has also conducted searches for binary inspiral sig- 
nals from primordial black holes (0.2-1.0 M ) in the halo 
of our galaxy [25| . plus more massive black hole systems 



where component masses are in the 3-20 M range [26[ . 

Bayesian inferential methods provide a means to use 
data from interferometric gravitational wave detectors 
in order to extract the parameters of a binary inspiral 
event. Markov chain Monte Carlo (MCMC) methods 
are a powerful computation technique for parameter es- 
timation within this framework; they are especially use- 
ful in applications where the number of parameters is 
large [27| . Nice descriptions of the positive aspects of a 
Bayesian analysis of scientific and astrophysical data are 
provided in [H [H [3(| ■ In previous work we have de- 
veloped MCMC routines for extracting five parameters 
associated with a binary inspiral event from data gen- 
erated by a single interferometric detector [3l|, [H [331 ] . 
Our MCMC code took data from a single interferometer, 
Fourier transformed it into the frequency domain, and 
then compared the result with 2.0 post-Newtonian (PN) 
stationary phase templates [34j . One of the new meth- 
ods that we implement in this current study, presented 
in this paper, is an MCMC routine that takes time do- 
main interferometric data, and compares it to time do- 
main templates that are 2.5 PN in phase, and 2.0 PN in 
amplitude; a trivial modification of the code (though not 
implemented in the study presented here) is to extend the 
accuracy of the signal waveforms to 3.5 PN in phase and 
2.5 PN in amplitude [H [H, H3, H, HI . A critical task 
for a world-wide gravity wave detection network will be to 
not only detect a binary inspiral signal, but to say where 
it came from. For this purpose the LSC has developed 
a coherent binary inspiral search pipeline [40l [4ll [42| . 
Coherent binary search pipelines and methods are also 
being developed within the Virgo collaboration [43| and 
TAMA [44[. Along similar lines, we have developed a co- 
herent MCMC parameter estimation routine, and in this 



Typeset by REVTeX 



2 



present paper we describe it and provide results from a 
test on simulated data. The simulations involve binary 
neutron star inspirals observed by three well-separated 
interferometers: the 4 km LIGO detectors at Hanford, 
WA and Livingston, LA, plus the 3 km Virgo detector 
in Cascina, Pisa, Italy. The synthetic data for the LIGO 
and Virgo detectors has Gaussian noise with power spec- 
tral densities (PSD) that match their target sensitivities 
[H, HI]. The MCMC code takes data from several in- 
terferometers, and estimates the two individual masses, 
time and phase at coalescence, distance to source, grav- 
ity wave polarization angle, angle of inclination of the 
binary system's orbital plane, and sky position in right 
ascension and declination. The additional parameters of 
polarization, inclination angle, and sky location can only 
be estimated accurately when data from more than one 
interferometer are considered; they also greatly inflate 
the parameter space and therefore complicate the analy- 
sis. MCMC methods have also been tested in a similar 
context to recover the nine parameters of a binary black 
hole coalescence in LISA data [13] ■ However, the problem 
setting is different (due to the different instrument, longer 
observation period, lower frequencies being investigated, 
and the 2.0 PN phase model), and MCMC techniques 
were employed rather for optimization than for integra- 
tion; the MCMC was also only applied on a subset of the 
parameters, while others were solved for analytically. 

The organisation of this paper is as follows. After a 
brief introduction to the analysis problem we describe 
our approach alongside more detail on the applied model 
in Sec. HB Sec. |TTT] provides practical directions how we 
implemented MCMC methods in order to analyze data 
within the described framework. Sec. IIVI eventually il- 
lustrates results of applications of our code to simulated 
data. We conclude the paper with a discussion and out- 
look in Sec. El 



II. ANALYSIS STRATEGY AND MODELLING 

A. Measuring gravitational waves 

In an inspiralling binary system, the two companions 
orbit around their centre of mass with decreasing or- 
bital distance and period, until the system eventually col- 
lapses. The gravitational radiation emitted by the system 
exhibits a 'chirp' form, that is, an oscillation of increas- 
ing frequency and amplitude. A laser interferometer is 
sensitive to space distortions along the directions of its 
two orthogonal arms, as it monitors the phase difference 
between the two laser beams that travel along the arms. 
A gravitational wave is a quadrupole wave that is char- 
acterized by its direction of travel, polarization angle, 
and its two polarization amplitudes. Its effect on a laser 
interferometer's measurements then is a linear combina- 
tion of the effects associated with the two polarizations, 
depending on the orientation of the interferometer with 
respect to the passing wave. Actual measurements, of 



course, are also exposed to noise. 

Measurements of a binary inspiral's chirp signal by a 
single interferometer will not be sufficient to infer all of 
the parameters that determine the signal's waveform and 
the interferometer's response. Measurements from sev- 
eral separate interferometers will, in general, be required 
to derive (for example) the wave's direction of travel by 
matching possible mutual orientations as well as different 
arrival times of the signals at the different sites. Combin- 
ing measurements from several interferometers will also 
enhance sensitivity and signal-to-noise ratio [40| . 

B. The Bayesian approach 

We apply a Bayesian approach to this inferential prob- 
lem, that is, the term 1 probability' 1 is used in a broader 
sense than in the more common ' frequentisV interpreta- 
tion [li, Hi, [4i| . Probability calculus here is applied to 
process and infer states of incomplete information that 
are reflected by probability distributions, and that are 
conditional on prior knowledge and/or the data at hand. 
This allows one to treat unknown parameters as random 
variables that follow a prior distribution representing the 
researcher's pre-experimental knowledge and uncertainty. 
The gain in information through observation of data then 
follows in a straight-forward fashion through the applica- 
tion of Bayes' theorem, yielding the posterior distribution 
of the parameters. The posterior distribution, which is 
essentially the product of the parameters' prior distribu- 
tion and the likelihood of the data, then poses the basis 
for inference [50 ]. 

Inference through the posterior distribution usually in- 
volves the solution of integrals, since one is typically in- 
terested in figures such as the marginal (posterior) ex- 
pectations of individual parameters, marginal (posterior) 
densities, or (posterior) probabilities of certain events, 
which are derived from the posterior distribution by in- 
tegrating over the parameter space. In many cases when 
analytic integration is not possible, numerical methods 
are employed, usually (pseudo-) stochastic techniques 
like Markov chain Monte Carlo (MCMC) methods that 
simulate random draws from the posterior ditribution, 
then allowing one t o appro ximate the desired integrals 
by sample statistics [23, [50] • 

C. Parameters 

The waveform that is measured at a certain interferom- 
eter depends on the characteristics of the inspiral event 
as well as the orientation of the source relative to the in- 
terferometer. The nine 'global' parameters determining 
the response of Earth-bound interferometers are: 

- individual masses (mi,m2 € 1R + ; mi < TO2), 

- luminosity distance (dz e K + ), 

- inclination angle (l g [0, 7r]), 
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- coalescence phase (0o S [0, 27r]), 

- coalescence time at geocenter (t c G It), 

- declination (8 G [—§,§]), 

- right ascension (a G [0, 2ir]) and 

- polarization (tp G [0, 7r]), 

the latter four of which affect the measurement at the 
I-th. detector in terms of the 'local' parameters 

- local coalescence time (t c G R) , 

- altitude (tfW G [0,tt]), 

- azimuth (ip^ 1 ^ G [0, 2%]) and 

- polarization (ip^ G [0, 7r]). 

These 'local' parameters are derived from the loca- 
tions/orientations of the source and the individual inter- 
ferometers with respect to each other. For more specific 
definitions and conventions see e.g. [5lj . In the follow- 
ing we will refer to the two parameter sets as the global 
parameter vector 

6»® = (mx,m 2 , d L , t, tf> , t c , 8, a, ip), (la) 

and the local parameter vector 

6& = (m u m 2> d L ,L,<t> ,t{ I) ,$ iI) ,<P iI) ,Tp (I) ) (lb) 

with respect to a specific interferometer /. Not all of 
the above parameters will usually be of primary interest; 
especially coalescence phase <fo, polarization ip or incli- 
nation l might be regarded as nuisance parameters. 

D. Network likelihood 

An interferometer's data output z is assumed to be the 
sum of the actual signal s(9), depending on the true pa- 
rameters 9, and (interferometer-specific) coloured noise. 
The (real- valued) data z and signal waveform s(9) en- 
ter the likelihood function in terms of their (complex- 
valued) Fourier transforms z and s(9), the noise is spec- 
ified through its power spectral density (PSD) S n . The 
likelihood function for a specific interferometer / then is 
(up to a normalising constant) proportional to the fol- 
lowing expression 

^>)o.e*p(-2f "Wy'V ) (2) 

[52| . For actual data, discretized and measured over a 
finite interval of length St , it is computed as the sum of 
squared and normalized differences between the Fourier 
transforms of the observed signal (z) and the signal tem- 
plate (s(9)) over the discrete set of Fourier frequencies 
{(£ X A/) : i L <i< iu}: 

data template 

oc exp ( ± V 

noise PSD 

(3) 



where ij x Aj and iu x Af are the lower and upper 
bounds of the examined frequency range. Note that, 
although not labeled as such here, data z, noise spec- 
trum S n , etc. are specific for the different interferome- 
ters /. Assuming that noise is independent across dif- 
ferent interferometers, the network likelihood then is the 
product of the individual interferometer likelihoods: 

C(9®) = ]JC iI} (9^). (4) 

E. Prior specification 

The prior information is specified in a straightforward 
fashion for the 'geometrical' parameters that determine 
the location and orientation of the inspiral event. A pri- 
ori, the event is assumed to be equally likely across all 
possible directions; this leads to a uniform prior for the 
right ascension a, and a prior density 

f(8) = ±cos(8) (5) 

that is proportional to the circumference of the corre- 
sponding 'circle of latitude' for the declination 8. Analo- 
gously, the prior density of the inclination i is 

/( i ) = ismW, (6) 

and the remaining parameters, polarization ip and 
coalescence phase <j)Q, have uniform priors. The prior 
specifications for these parameters may also be regarded 
as Maximum Entropy choices [49l [531 ] . 

The coalescence time t c is assumed to be known in 
advance up to a certain accuracy through preprocessing 
of the data [13, H3, ; for now we set its prior to be 
uniform across ±5ms around the true value, which of 
course is known for simulated data. 

The joint prior of the remaining parameters, masses 
mi, m 2 and luminosity distance d^, is set in order to re- 
flect the distribution of parameters given the event has 
been detected in the first place. Initially, the prior for the 
two inspiral companions' individual masses is assumed 
to be uniform between 0.6 and 3.0 M (solar masses: 
M w 2x 10 30 kg), which effectively covers the range of 
values expected for neutron stars. The prior for rf^ is 
derived from the assumption that inspirals happen uni- 
formly across space, so that P(dL < i) « i 3 . So far, this 
leads to an improper distribution (that has an infinite 
integral) . But inspiral signals obviously cannot originate 
from arbitrarily great distances, since at some point their 
signals become too faint to be detected. We incorporated 
this restriction by taking into consideration the detection 
probability D, which we assume to depend on the signal's 
amplitude, which is roughly proportional to 

, y , , . / Jm\m2 \ 
A(m 1 ,m 2 ,d L ) = in — -j- , (7) 

\d L (mi + to 2 ) 6 / 
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neglecting for simplicity the effects of orientation param- 
eters (.4 is actually the logarithmic amplitude) [33|. We 
could have set a threshold amplitude below which in- 
spirals would be assumed undetectable, but favoured a 
'smoother' transition that does not explicitly apply zero 
probability to parts of the parameter space. Instead we 
model the dependence between signal amplitude and de- 
tection probability using a (sigmoidal) logistic function 
of the form 



D a ,b{x) 



1 



l + exp(^) 



(8) 



whose parameters a and b are set so that D a j,{xu) = 1 — p 
and D a ,b(xL) = p, for some upper and lower reference 
points xjj and Xl, and some < p < 0.5 (e.g. p := 
0.1). So Xjj denotes the amplitude at which the detection 
probability reaches 1 — p, and Xl is the amplitude at 
which the probability falls below p. In order to fit D 
through these points, its parameters are set to 



xl + x v 



and 



xy - x L 
21nW 



(9) 



So the density of the resulting (proper) prior distribution 
of individual masses and distance is 

f(m 1 ,m 2 ,d L ) = cx d 2 L x D a . b (A(mi,m 2 ,d L )) (10) 

for some normalizing constant c € R + [33] . For the exam- 
ples shown here, we set xu ■= .4(2. 0M o , 2.0M o , 45Mpc), 
x L := .4(2.OM ,2.OM ,5OMpc) and p := 0.1, so an 
optimally oriented 2.0-2.0 M inspiral is assumed to be 
detectable out to 45 and 50 Mpc with probabilities of 
90% and 10%, respectively. 

As a 'side effect' of this prior definition, larger masses 
have a greater prior probability, since inspirals involving 
large masses may originate from farther distances while 
low-mass inspirals need to be close in order to be observ- 
able at all. This feature is also known as the Malmquist 
effect; incorporating it into the prior will compensate 
for selection bias that would otherwise affect the results 
[H, H3- The definition of priors, especially for coales- 
cence time t c , individual masses mi and 771,2 and their 
relation to the luminosity distance d^, and possibly also 
for the sky location (6, a), may be refined at a later 
stage when e.g. some substantiated knowledge is avail- 
able about the performance of the upstream detection 
pipeline, which might provide rough estimates of some of 
the parameters together with the detection [12, [H, [55| . 
For now we aim for simple and general formulations. 



III. IMPLEMENTATION 
A. General 

In order to analyze data in terms of the above frame- 
work, we implemented an MCMC sampler in C that is 



supposed to both find the global mode(s) of the posterior 
distribution and then 'explore' the distribution, i.e. simu- 
late random draws from the posterior. Data is imported 
from the Frame format using the Frame Library [58| . 
Prior to the analysis, the data is filtered and downsam- 
pled by a factor of 4 [5i|[6(|. (Pseudo-) random number 
generation within the MCMC sampler was implemented 
using Randlib [6l|. 

The MCMC sampler writes only every 25th of the 
drawn samples to a text file, in order to reduce the ef- 
fects of subsequently correlated samples, and also to keep 
the data volume at a reasonable level. The MCMC out- 
put then is imported into R, a statistical software, for 
eventual analysis [62| . 

The marginal densities that are shown in this paper 
are kernel density estimates 63]. Two-dimensional den- 
sities are illustrated by greyscale plots (with darker areas 
corresponding to greater densities), and in addition the 
contour line enclosing the most probable region (accumu- 
lating 95% probability) is shown. 



B. Likelihood implementation 

In order to compute the coherent network-likelihood, 
first the individual interferometer-likelihoods need to be 
determined. The primary 'ingredients' for the interfe- 
rometer-likelihood are 

- the Fourier-transform of the data z, 

- the noise spectrum S n , 

- the 'local' parameter set 9^ and 

- the (Fourier-transformed) signal template s(9^), 

the first two of which only need to be determined once at 
the very beginning of the analysis, while the latter two 
(in general) need to be re-computed for each likelihood 
evaluation. 

For all (discrete) Fourier-transformations we use the 
freely available FFTW library [64[. The noise spectrum 
is estimated from a section of data that is disjoint from 
the actually analyzed data set 65]. In order to mini- 
mize undesirable leakage effects, the data is 'windowed' 
before Fourier-transformation; using a Hann window for 
spectrum estimation, and a Tukey window for data and 
template transformations [66| , 

Internally, along with the noise spectra, data Fourier 
transforms etc. corresponding to each of the interfer- 
ometers /, a set of vectors defining the interferometer's 
location and orientation is stored. This allows to de- 
rive the interferometer-specific parameters (local coales- 
cence time t^P, altitude i9^, azimuth ip^ and polariza- 
tion ip^) with respect to the galactic and Earth coor- 
dinate systems from the global parameter vector 9 9 via 
vector operations like rotations, orthogonal projections, 

etc. M,m H. 
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C. Time-domain waveform generation 

Template waveforms s(9) are generated in the time- 
domain and then (numerically) Fourier transformed to 
the frequency domain. Here we used waveform approxi- 
mations that are 2.5 PN in phase and 2.0 PN in ampli- 
tude. The rather complex expressions for these are omit- 
ted here and can be found in [HT(] . We preferred work- 
ing in the time-domain, since frequency-domain tem- 
plates might introduce discrepancies because they are ex- 
act analytic Fourier transforms of the 'parametric' wave- 
forms, while the actual data is of finite extent and af- 
fected by leakage introduced through the numerical dis- 
crete Fourier transformation. When using time-domain 
templates and Fourier transfoming these, the resulting 
frequency-domain templates are the more accurate match 
to the Fourier-transformed data. Another advantage of 
generating templates in the time-domain is the availabil- 
ity of a wider range of signal waveforms; the extension to 
higher PN appoximations ( e.g. 3.5 PN phase and 2.5 PN 
amplitude [H, [H, H3, [H, [39|) or the consideration of 
additional parameters (e.g. spin effects 70]) would be 
straightforward to implement. 



Parallel tempering is a special case of the 'Metropolis- 
coupled MCMC (MCMCMC) algorithm, in which sev- 
eral 'tempered' chains are run in parallel, each having a 
different (constant) temperature. So the tempering does 
not vary over time, but instead is realised across paral- 
lel chains, with additional proposals allowing for swap- 
ping between chains. Instead of sampling from the regu- 
lar posterior distribution with density function / (which 
is essentially the product of prior it and likelihood C: 
f(0) oc n(9)C(9)), the tempered chains sample from a 
modified distribution 



f T (9) ex tt(9)C(9)t, 



(12) 



where T > 1 denotes the 'temperature', and for 
which in the extreme cases T — 1 and T — > oo the 
tempered distribution fx equals posterior and prior 
respectively. Chains running at higher temperatures can 
be considered as sampling from a 'relaxed' or 'widened' 
posterior which is then used as proposal distribution 
(through the swapping between chains) for 'cooler' 
chains, thus improving both convergence and mixing. 
The draws from the 'coolest' chain with T = 1 are the 
only ones that are eventually used for posterior inference. 



D. MCMC implementation 

In order to enhance the MCMC sampler's performance 
we applied reparametrisations to some of the parame- 
ters. The individual masses (mi, mo) are highly cor- 
related in their posterior distribution [181 ]. making sam- 
pling from the original parameters extremely difficult. 
Re-expressing the masses in terms of chirp mass m c and 
the (symmetric) mass ratio r\, where 



m c = 7 \T7K and V = 7 v7 

{mi+m 2 ) L > b {mi+m 2 y 



(11) 



yields a posterior that is much easier to sample from. We 
then reparameterized the luminosity distance from to 
ln(c?i), which implicitly yields an unbounded parameter 
space and proposal step widths that are proportional to 
the current distance oIl- Declination 8 and inclination i 
were transformed to sin(<5) and cos(i), which leads to 
uniform prior distributions over the new parameters. 

The MCMC algorithm was implemented as a Metrop- 
olis-sampler [13, H3] that was extended to a parallel tem- 
pering algorithm. The idea of tempering is to consider 
a smoothed ('tempered') version of the actual objec- 
tive function (here the posterior distribution), or, anal- 
ogously, do its exploration (optimization, MCMC sam- 
pling,...) following 'relaxed' rules. In optimization con- 
texts, the tempering is often faded out over time, the 
result being a 'simulated annealing' algorithm, which 
starts off at a high temperature in order to find the 
global mode amongst other minor modes, but eventu- 
ally ends up optimizing the original objective function. 



The starting values for the sampler are determined us- 
ing importance resampling (Hoj . In a simplified setting of 
the problem (considering only 5 parameters and one in- 
terferometer [33j), this was sufficient to yield reasonable 
posterior samples that were close enough to the main 
mode so the sampler then converged reliably and fast. 
Due to the much larger parameter space and the com- 
putationally more expensive likelihood, ensuring conver- 
gence through good starting values is not feasible any 
more. Instead, convergence is now supported through 
the use of parallel tempering, which enables the sampler 
to cross gaps between (local) modes and eventually find 
the posterior's global mode(s). 

As proposal distribution for the MCMC sampler we 
used a multivariate Student's i-distribution with 3 de- 
grees of freedom. In addition to these 'regular' proposals, 
sometimes draws from the prior are proposed for some 
parameters in order to enhance convergence, or steps to 
'related' parts of the parameter space, like a step from 
phase 0o to </>o ± 7T, which corresponds to an (almost) 
equally likely parameter combination if the two masses 
are (almost) equal. Proposals like these are valid as long 
as a certain symmetry is maintained, i.e. every proposed 
step is as likely as the reverse step; otherwise one would 
need to switch to a Metropolis-Hastings sampler that is 
able to account for asymmetric steps [27J, |50j . 



E. Signal-to-noise ratios 

The signal-to-noise ratios (SNR) stated in subsequent 
sections are defined as follows. The interferometer- 
specific SNR of a certain signal s(9®) received at interfer- 
ometer / and embedded in noise with spectral density S n 
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is defined as: 

[32| . We computed it — in analogy to Eq. ([3]) — over the 
same frequency range that is relevant for the likelihood. 
The network SNR then is defined as: 




IV. EXAMPLE APPLICATION USING 
SIMULATED DATA 

A. Overview 

In the following two sections we will illustrate results 
generated by our MCMC implementation on simulated 
data. Firstly, in section IV. B the evidence on parame- 
ters in the measured data is illustrated for one simulated 
signal by looking at the results of an MCMC run in detail. 
In the following section IV. C the effect of varying signal 
properties on the evidence in the data is demonstrated 
and some peculiarities are pointed out. The posterior dis- 
tributions of parameters are compared across simulated 
signals that are observed at different distances (and with 
that, SNRs), but which are otherwise identical. While 
the individual SNRs (at each interferometer) for each of 
these examples are similar, a more extreme example with 
an almost zero SNR at one of the interferometers is con- 
sidered as well. 



B. Recovering the inspiral's parameters 

We simulated an inspiral event that is measured at 
three interferometer sites, namely Hanford (LHO), Liv- 
ingston (LLO) and Pisa (V). Due to their different noise 
characteristics, the frequency ranges of the likelihoods 
were set to 30-1600 Hz for the Virgo observatory, and 
40-1600 Hz for the other two LIGO interferometers (cp. 
Eq. The amount of data to consider was set to 

be the 40 seconds before coalescence for Virgo, and 20 s 
for the other two. This matches the time an inspiral of 
this kind spends emitting radiation within the above fre- 
quency ranges, and would in a realistic search need to 
be set either with respect to worst-case considerations, 
or based on prior information supplied by the detec- 
tion pipeline. The original sampling rates of the data 
were 20 000 Hz (V) and 16 384 Hz (LHO, LLO), and the 
signals were superimposed with Gaussian noise match- 
ing the corresponding interferometer's design sensitivities 
[45I |46| . The example inspiral had parameter values of 
c?l = 10 Mpc for the distance and masses of mi =1.5 M 



and TO2 = 2.0 M (chirp mass m c — 1.505 M , mass 
ratio 77 = 0.245). The resulting SNRs at the three inter- 
ferometer sites are: LHO: 16.4, LLO: 21.2 and V: 12.6 
(network SNR: 29.6). 

Six parallel MCMC chains were run within the par- 
allel tempering scheme. With this amount of data the 
MCMC code generated some 80 samples per minute on a 
3.2 GHz Pentium desktop PC, so considering the parallel 
chains (six) and the thinned output (every 25th), an ac- 
tual posterior sample is generated every 113 seconds. The 
first chain of the parallel tempering MCMC sampler con- 
verged after some 75 000 iterations, after 'thinning out' 
of the samples and discarding the burn-in phase, the re- 
sulting posterior sample size was 12 500 samples. 

FIG [1] shows marginal posterior densities of the nine 
parameters for our example problem, and Table U lists 
some numerical posterior estimates. Although correla- 
tions between parameters are already greatly reduced 
through the reparametrisation, some correlation still re- 
mains. FIG. illustrates correlations between two such 
pairs of parameters, one can see that the 'new' mass pa- 
rameters m c and r\ are still dependent (though orders of 
magnitude less than m\ and mi were), and also that the 
uncertainty in the luminosity distance d^ is tied to the 
uncertainty in the inclination angle 1, since these two pa- 
rameters can compensate or mimic each other's effect to 
a certain degree. 

Distributions of other variables derived from the pa- 
rameters could be investigated, their distributions then 
dependin g on the joint distribution of the involved pa- 
rameters [33J. Examples would be the individual masses 
(mi, TO2), or the total mass m t = mi + WI2; the distri- 
bution of mt is narrower than those of both mi and m,2, 
due to their negative correlation (cp. Table [I]). 



TABLE I: Some key figures summarizing the marginal pos- 
terior distributions of individual parameters, where meaning- 
ful. Mean and median characterize the distributions' cen- 
ters. Given the observed (simulated) data, the parameters 
fall within the central posterior intervals with 95% probabil- 
ity. The true parameter values used to generate the data are 
shown as well. 





mean 


median 


95% c.p.i. 


true 


unit 


m c 


1.5044 


1.5044 


[1.5039, 1.5048] 


1.5047 


M 


V 


0.2418 


0.2417 


[0.2380, 0.2460] 


0.2449 




tc 


12.3445 


12.3445 


[12.3440, 12.3450] 


12.3450 


s 


d L 


8.89 


8.29 


[6.25, 13.1] 


10.00 


Mpc 


6 


-29.76° 


-29.77° 


[-30.74°, -28.84°] 


-29.00° 




a 


17 h 45.0' 


17 h 44.9' 


[17 h 42.1', 17 h 48.9'] 


17 h 45.0' 




L 


0.911 


0.995 


[0.194, 1.354] 


0.700 


rad 


mi 


1.446 


1.442 


[1.389, 1.526] 


1.5 


M 


rn 2 


2.080 


2.084 


[1.965, 2.170] 


2.0 


M 


m t 


3.526 


3.527 


[3.490, 3.559] 


3.5 


M 
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FIG. 1: Marginal posterior densities of the inspiral's parameters for our example problem. Dashed lines indicate the true values. 




1.504 1.505 o t/4 it/2 

chirp mass m c (M Bun ) inclination i (rad) 



FIG. 2: Marginal joint posterior distributions of two pairs of 
parameters. Dashed lines indicate the true values. 



C. Results for varying signal characteristics 

We performed additional runs with varying settings of 
the 'true' parameters of the simulated signal. As one 
would expect, the precision of parameter estimation is 
proportional to the signal's strength; Table [IT] shows the 
standard deviations of some of the parameters' posterior 
distributions. The posterior is narrowest for a close- by 
inspiral of high masses, and gets wider for both lower 
mass or greater distance. 

These results are in agreement with earlier estimates 
of the accuracy to be expected from such parameter esti- 
mates [ll]. The great difference in relative accuracies of 
parameters related to phase evolution (like chirp mass m c 
and reduced mass u = mi , m2 = m t ri) versus those affect- 
ing the signal's amplitude (like distance d^) is confirmed, 
and the correlation between m c and fi is verified as well. 

At decreasing SNRs, certain parameters cannot be de- 



termined unambiguously any more. One example is the 
inclination angle which still has a 'well-behaving' pos- 
terior distribution at 10 Mpc distance (see FIG. For a 
weaker signal originating from 30 Mpc distance, the dis- 
tribution then turns bimodal (FIG. [3]). The 'orientation' 
of the inclination angle is not clear any more, the result 
being two roughly equally likely 'mirror image' solutions 
with P(t < f ) w | w P(t > f ). Note that the two solu- 
tions l and 7r — l correspond to opposite orbital directions 
(clockwise/counterclockwise), as seen from Earth, which 
might be of minor interest anyway. 

The sky location's posterior also exhibits multiple 
modes for this weaker signal (FIG. [3]). This illustrates 
some potential pitfalls of Maximum-Likelihood (ML) or 
Maximum-a-Posteriori (MAP) methods; these would ad- 
vise picking the highest of the several modes, which might 




it/2 „ 18 h 17 h 16 h 

inclination i (rad) right ascension a 



FIG. 3: At greater distance the 'orientation' of the inclina- 
tion angle t cannot be resolved any more, both directions are 
roughly equally likely (P(t < f ) « ± « P(t > f )). At the 
same time, with the lower SNR the sky location's posterior 
turns multimodal. (Dashed lines indicate the true values.) 
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TABLE II: Individual and total SNRs for different signals, and some characteristics of the resulting posterior distributions. The 
accuracy of some of the parameters is illustrated by the posterior standard deviations for (5, a), t c , d^, m c and fi (percentages 
refer to the true value). The correlation coefficient for m c and /i shows the (posterior) interdependence between the two 
parameters. Our results are consistent with those presented in |l8| ]. 



masses distance signal-to-noise ratios posterior standard deviations 



mi-m,2 


d L 


LHO 


LLO 


V 


network 


(6,af 


tc 


d L 


m c 




Cor(m c ,^) 


1.5-2.0 M 


10 Mpc 


16.4 


21.2 


12.6 


29.6 


0.011 rad 


0.26 ms 


20% 


0.016 % 


0.35 % 


0.95 


1.5-2.0 M 


20 Mpc 


8.2 


10.6 


6.3 


14.8 


0.030 rad 


0.49 ms 


25% 


0.031 % 


0.69 % 


0.94 


1.5-2.0 M 


30 Mpc 


5.5 


7.1 


4.2 


9.9 


0.207 rad 


1.04 ms 


25% 


0.074 % 


1.33 % 


0.91 


2.0-2.0 M 


10 Mpc 


18.4 


23.9 


14.1 


33.3 


0.008 rad 


0.14 ms 


14% 


0.009 % 


0.14 % 


0.80 


2.0-2.0 M 


20 Mpc 


9.2 


11.9 


7.1 


16.7 


0.017 rad 


0.28 ms 


18% 


0.014 % 


0.23 % 


0.73 


2.0-2.0 M Q 


30 Mpc 


6.1 


8.0 


4.7 


11.1 


0.026 rad 


0.42 ms 


21% 


0.021 % 


0.37 % 


0.78 



a spherical standard deviation [7ll 



just be the narrowest one, but not necessarily the most 
likely. If one then proceeded by extrapolating the cur- 
vature at that mode and deriving error bounds from the 
Fisher Information matrix, the resulting estimates might 
not only be far off, but also associated with overestimated 
accuracies. 

We also tried MCMC runs with a modified prior set- 
ting; we extended the prior for the coalescence time t c 
from its original range of ±5ms around the true value to 
±27ms, allowing for an additional margin of 22ms, which 
is the time it takes a gravitational wave to travel from 
Earth's surface to its center. This setting reflects the 
case where the inspiral detection pipeline received trig- 
gers from less than three interferometer sites, so the sig- 
nal's arrival time at the geocenter could not be estimated 
to greater accuracy in advance. The MCMC algorithm is 
still able to find the mode in the enlarged time parameter 
range, but takes more iterations to converge. 

One scenario in which such an approach would be 
necessary is when the SNR for one of the interferom- 
eters is almost zero. In such a case the data from 
the interferometer under consideration also would not 
(directly) contribute to the estimation of phase- and 
frequency- related parameters, but would still carry in- 
formation about amplitude-related parameters — by im- 
plicitly 'ruling out' those parameter combinations that 
would have resulted in a response at that interferome- 
ter. FIG. H] shows the sky location posteriors for such a 
signal, a 1.5-2.0 M Q inspiral at 10 Mpc distance, where 
the SNRs at the three interferometer sites are: LHO: 9.6, 
LLO: 13.9, V: 0.18 (network: 16.9). Including the data 
from the third interferometer (with almost zero SNR) 
into the analysis still yields a much more accurate esti- 
mate of the sky location. Table ITTT1 compares the resulting 
parameter accuracies of these two settings. The posteri- 
ors for sky location (<5, a) and coalescence time t c , which 
are closely related, are much narrower when the Virgo 
data is considered in the analysis, while estimates of the 
rather phase- and frequency-related parameters m c and fi 
do not gain from the additional information. 

On the one hand, not only a high (total) SNR is desir- 
able but also one that is rather 'evenly spread' over differ- 
ent interferometers. On the other hand even a near-zero 




19 h 18 h 17 h 19 h 18 h 17 h 

right ascension a right ascension a 

FIG. 4: Even if the SNR at one of the interferometers is 
almost zero, it still contributes to estimates' accuracies — the 
posterior is much narrower if its data is included (left plot) 
than if it is omitted (right plot). (Dashed lines indicate the 
true values.) 



SNR at one of the interferometers does not make its mea- 
surement useless. Inference on different parameters will 
be affected to different degrees by such an unbalanced 
SNR arrangement. 

TABLE III: Relative accuracies of different parameters (in 
analogy to Table [TTJ) when considering / not considering the 
Virgo data (where the example situation is such that the SNR 
is nearly zero). 



Virgo data... 


(S, a) 


tc 


dL 


m c 




A 4 




...included 


0.071 rad 


0.81 ms 


21 % 


0.023 


% 


0.33 


% 


...excluded 


0.150 rad 


2.38 ms 


23% 


0.022 


% 


0.31 


% 



V. DISCUSSION 

We have developed a new MCMC program for esti- 
mating the nine parameters associated with an inspiral 
of compact binary objects from the data coming from a 
network of gravitational wave interferometers. The de- 
termination of the sky location of the source is an im- 
portant consequence of the procedure. Numerous new 
features have been implemented in this binary inspiral 
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MCMC. The MCMC uses waveform approximations that 
are 2.5 PN in phase and 2.0 PN in amplitude [5l| (and by 
the time of final submission of this paper a version using 
waveforms that are 3.5 PN in phase and 2.5 PN in ampli- 
tude [H, [H, [13, [H, HI] was running as well). The data 
from multiple interferometers (two or more) are coher- 
ently analyzed in order to produce posterior probability 
distributions for all nine parameters. 

Advanced MCMC techniques were implemented in our 
program in order to maximize the efficiency of con- 
verging to the correct parameter values in the large, 
9-dimcnsional, parameter space. The initial parameter 
values for the sampler were determined using importance 
resampling 50]. We recently extended (though not with 
the results presented in this paper) the parallel temper- 
ing algorithm to an Evolutionary MCMC algorithm [72] • 
This MCMC variety implements proposals that are moti- 
vated by genetic algorithms [73], and so recombinations 
of parameter samples from different MCMC chains are 
used as proposals in order to improve convergence and 
mixing. 

Another current related research effort is the applica- 
tion of a version of this MCMC code to burst waveforms. 
This problem is by orders of magnitude computation- 
ally less expensive, due to the much shorter duration 
of the signals. But it appears that on the other hand 
convergence, i.e. finding the main posterior mode in the 
parameter space, still poses a major problem. The theo- 
retical background of the various potential burst sources 
is rather vague, so realistic waveforms and sensible spec- 
ifications of parameterisations and priors also need to be 



identified. The results of this study on the MCMC pa- 
rameter estimation of burst signals using the coherent 
analysis of multi-interferometer data will be presented in 
a forthcoming publication. 

We are also extending the MCMC techniques used in 
this study to the application of data analysis for LISA 
detection of binary inspiral signals. While our present 
program coherently analyzes data from multiple ground 
based interferometers, we have found it is a straightfor- 
ward extension of the code so that we can coherently 
analyze the time delay interferometry data from LISA. 
These results are also forthcoming. 

Presently LIGO is at its target sensitivity. Virgo is 
fast approaching its design sensitivity. Using the LIGO- 
Virgo network it will be possible to observe neutron star 
binary inspirals out to a distance of 35 Mpc [45|, |4(| . A 
detection of such an inspiral could occur at any time 0] • 
As displayed in this paper, our MCMC routine is capa- 
ble of coherently analyzing the data from the multiple 
interferometers, and then using it to estimate the nine 
parameters associated with such a signal. 
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